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1. Introduction 

A quantitative first-principles theory of atoms, molecules, and solids re- 
quires an accurate description of the electron-electron correlation. Fahy's 
[1], Mitas' [2], and Umrigar's [3] contributions to this book suggest that 
Variational or Diffusion Monte Carlo simulations are rapidly approach- 
ing the required physical or chemical accuracy. However, at least at their 
present pioneering stage, such simulations are still limited to relatively sim- 
ple systems, whose correlation energy is a small fraction of the total energy, 
and, more significantly, whose main physics and chemistry can be qualita- 
tively understood without any reference to many-body theory. By contrast 
(and by definition), strongly correlated electron systems are those for which 
even a qualitative understanding of the main physics and chemistry is be- 
yond the independent-electron picture; in this case, as also pointed out by 
Koch in his contribution to this book [4], rather crude models of real in- 
teracting electrons, like the Hubbard model, can be on one hand the only 
feasible target, and, on the other, a clever way to capture the essential 
physics of a whole class of materials without going into too many realis- 
tic details [5]. However, apart from very special limiting cases, even simple 
models of interacting electrons cannot be solved exactly, because, in spite 
of many severe approximations, they still give rise to an astronomically 
large Hilbert space; that's why quantum Monte Carlo methods, based on 
the idea of stochastic sampling, remain a natural option. In this new con- 
text, however, the absence of a preliminary qualitative understanding of the 
electronic correlations makes the choice of the trial wavefunction (Sect. 2), 
and possibly even of the size of the simulation box (Sect. 4), a much more 
crucial and less straightforward part of the investigation. A recent work by 
Hood et al. [9] gives the occasion of an instructive glance at the state of 
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the art, illustrating the different theoretical challenges involved in today's 
QMC studies of "realistic" and "model" systems. Hood's variational Monte 
Carlo investigation of bulk silicon, based on Jastrow-Slater wavefunctions 
and a 54-atom 3D lattice, is more than adequate to obtain an essentially 
exact ground state, thereby explaining the great success of approximate 
local-density functionals for this kind of materials (similar considerations 
apply e.g. to s-p bonded materials, like boron nitride [1, 10]). By contrast, 
when investigating the Hubbard model, the extensive Fixed-Node Monte 
Carlo simulations by Cosentini et al. [11], based on the nodal structure 
of a correlated Gutzwiller-Slater wavefunction and a (proportionally much 
larger) 256-atom 2D lattice, may in fact approach the true ground state of 
the infinite 2D Hubbard lattice with a much better accuracy than most pre- 
vious attempts, and yet may not be sufficient to proclaim the ultimate word 
on phase separation. Rather than giving a final answer to this really tough 
problem, the study presented in the second part of this chapter (Sects. 3 
and 4) thus serves the more modest purpose of illustrating a particularly 
challenging application of the Fixed-Node quantum Monte Carlo method, 
recalled in the next Sect. 2. 

2. Fixed-Node QMC for lattice fermions 

In what follows we briefly illustrate some aspects of the Fixed-Node QMC 
scheme actually followed in the study of the one-band Hubbard model of 
the subsequent Sects. 3 and 4. We omit all demonstrations, which can be 
found either in other chapters of this book or in the original papers (in 
particular Ref. [12] for the power method and Ref. [13] for the fixed-node 
approximation) . 

2.1. POWER METHOD ON A LATTICE 

Lattice models of quantum many-body systems (see e.g. Ref. [14]) are nor- 
mally based on the preliminary choice of a tiny single-body basis set at- 
tached to each site of a lattice (for example, just one orbital per site in the 
simplest Hubbard model, or one spin per site in the Heisenberg model). On 
a finite lattice of N s sites, the corresponding many-body Hilbert space is 
also discrete and finite, but its size is a very rapidly growing function of N s 
(for example, it's exponential for Heisenberg, combinatorial for Hubbard). 
Since the purpose of lattice hamiltonians is to model infinite solids, small- 
N s lattices (which can be either analytically solved or exactly diagonalized 
on the computer) may be interesting as a benchmark for new theories and 
numerical methods (see e.g. Table 1), but the real challenge is to reliably 
extrapolate the behavior of the infinite model solid from lattices of reason- 
ably large N s . As mentioned, this amounts to dealing with finite but huge 
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Hilbert spaces, where quantum Monte Carlo methods are often the only 
available alternative. 

2.1.1. A sparse matrix whose elements are all positive 

In this context it is possible to project the ground state out of an arbitrary 
initial state (not orthogonal to it) by the "power method" - the repeated 
application of some operator Q simply related [15, 16] to the hamiltonian 
Ti (often conventionally called Green's function) - provided that, in the 
(huge) many-body Hilbert space of the lattice model, one can choose a rep- 
resentation such that (a) each of the basis states \X) is connected only to a 
few other basis states \X') by the operator Q (or, in other words, most off- 
diagonal matrix elements of Q are zero and the matrix {X'\Q\X) is sparse), 
and (b) all the nonzero matrix elements of Q are positive (X'\Q\X) > 0. 
As described by Nightingale's contribution [17], the condition (a) is needed 
for the action of the hamiltonian TL (or related [15] operators Q) on some 
arbitrary state to be computationally feasible - the number of operations 
involved must be much smaller than the size of the Hilbert space [18]; the 
condition (b) implies that, in the chosen representation, (X'\Q\X) can be 
split into a Markov matrix and a positive weight, and stochastic approaches 
become possible (see e.g. Ref. [17]); it also implies that the ground-state 
wavefunction [16], in that particular representation, is positive everywhere 
(ATl^o) > 0, which rules out systems with three or more fermions (sign 
problem). In other words, the power method presented here does not apply 
to the ground state of lattice fermions unless some approximation (recalled 
in Subs. 2.2) is adopted; it instead directly applies to those lattice mod- 
els (like e.g. the Heiseberg model on a square lattice [12]) for which both 
conditions (a,b) are satisfied. 

2.1.2. Markov chain (no weights) 

Let us first suppose for purely pedagogical reasons that, besides (a) being 
sparse and (b) having all elements greater or equal to zero, the relevant 
matrix also (c) has the additional, remarkable property ^x'(X'\G\X) = 1. 
If (c) holds, then (AT'lC/jAr) can be directly identified with a probability 
(Markov) matrix Px'X, whose maximum right eigenvector |0 max ) (i-e. the 
right eigenvector corresponding to its largest eigenvalue A max = 1) can be 
first approached, and then sampled, by a (sufficiently long) single random 
walk, which can be generated on the computer by e.g. the following simple 
rule: 

1) take an initial state |X ) 

2) find all the N Q states \X') for which P x , Xo = (X'\g\X ) ^ (including 
X'=X Q ) 
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3) evaluate Px'X a = (X'\Q\X ) > for each of the above states, obtaining 
a list of N probabilities which add up to 1 

4) align them as N a adjacent segments adding up to form a unit segment 

5) take a random number between and 1 and, depending where it falls 
on that unit segment, choose the corresponding state \X') as the new 
initial state \X\) 

6) go back to the first item of this rule and iterate K times this process 
(with very large K) 

The above rule amounts to starting from an initial state \X Q ) and to re- 
peatedly (K times) multiplying from the left the matrix Px'x- The matrix- 
vector multiplication is implemented each time in a statistical sense: for 
example in the first step, rather than yielding a linear combination of all 
the N Q states \X') with coefficients Px'X , the statistical result of the 
multiplication is just one of them, chosen at random among all the pos- 
sible N Q states \X') according to the probabilities Px'X - After a good 
number n of such stochastic matrix multiplications, the random sequence 
\X n+ i), \X n+ 2), ■ ■ ■ \Xr) starts to be distributed according to \4> ma , x ). If the 
condition (c) holds, then the repeated multiplication of Px'x and the re- 
peated application of the operator Q are the same thing, \4> max ) is pro- 
portional to l^o), and the random walk \X n+ i),\X n+ 2), ■ ■ ■ \Xk) can be 
directly used for statistical averages over the ground state. 

2.1.3. Single random walk (weights) 

Unfortunately general Markov matrices Px'x are not symmetric (that's 
why one has to distinguish between left and right eigenvectors), and H or 
related hermitian operators Q do not satisfy the condition (c) of previous 
Subs. 2.1.2: one has, instead, J2x'(X'\G\X) = bx ^ 1. However, if the 
property (b) does hold, then bx > 0, and we can still decompose (X'\Q\X) 
into the product of a normalized probability matrix Px'x an d a positive 
weight bx- Then, by the simple rule illustrated above (steps 1-6) we can 
still generate a random walk which samples the maximum right eigenvector 
°f Px'x = (X'\Q\X)/bx, except that now the eigenvector |0 max ) is no 
longer proportional to |vP )j t ne ground-state eigenvector of 7i. However, 
if we keep record not only of the state \Xj), but also of the weight bj = 
bxj = J2x'(X'\Q\Xj) (which must be calculated at each step because it's 
needed to obtain a Markov matrix P from the original matrix Q), then 
the random walk (n<i<K) which samples \4> max ) can a ls° be used 

to sample £/ L |</> m ax): one only needs to associate an L-th order cumulative 
weight Wi(L) = bi_Lbi_L + ibi-L+2 ■ ■ ■ h-i to each state |Aj) of the sequence, 
and then consistently include these weights in the averages [12]. Since the 
repeated application of Q tends to project the ground-state component out 
of l^max), a sufficiently large L could allow ground-state averages from a 
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single random walk [19]. Unfortunately this is not possible because, when 
L is large enough for Q L to project the ground state out of |</>max), then the 
weights Wi(L) - as suggested by their product form and argued in better 
detail elsewhere [12] - undergo wild fluctuations, which grow exponentially 
with L. So statistical averages are accompanied by a diverging variance and 
become pointless. 

2.1.4. Statistical averages, many walks, reconfiguration 
The standard solution to the problem of a diverging variance, which occurs 
with a single random walk, is to propagate many simultaneous random 
walks. These walks propagate independently most of the time, but not all 
the time (else they would be equivalent to a single, long walk): to prevent 
individual weights from getting too small or too large and the variance 
from blowing up, the different walks must periodically undergo some sort 
of global reshoveling of states and weights. As recently emphasized by Ca- 
landra and Sorella [12], it's this reconfiguration process which makes the 
difference. The results of Sects. 3 and 4 are actually based on their new 
reconfiguration scheme. A thorough discussion of the underlying theory, 
and even a complete presentation of their recipe, is beyond the scope of 
this chapter; it can be found by the interested reader in their recent, origi- 
nal paper [12]. But our numerical experience may be worth mentioning: in 
practice the Calandra-Sorella reconfiguration scheme allowed us to repro- 
duce with a relatively small fixed number (100 ~ 200) of walkers equally 
accurate results as those obtained by means of standard branching schemes 
[20] with more than 2000 walkers. 

2.2. THE SIGN PROBLEM 

As mentioned in the previous Subs. 2.1, the power method can work as a 
ground-state [16] projector only if all the matrix elements of the operator 
Q are positive [15], or, equivalently, if the off-diagonal matrix element of 
the hamiltonian are all negative [15, 16]. 

2.2.1. Fermion nodes and other possible sources 

The above condition on matrix elements can be violated for various reasons: 
(i) particle statistics: the ground-state wavefunction of identical fermions 
changes sign under their exchange, and must be negative somewhere; (ii) 
one-particle effects: in certain lattice models, some hopping terms have 
the "wrong" sign (on continuum "realistic" models a similar problem may 
arise from nonlocal pseudopotentials [21] as illustrated by Mitas in this 
book [2]); (iii) other effects, like the frustration of a quantum spin system 
on a triangular lattice [22] ; (iv) combinations of the above. From this point 
of view our simplest 2D one-band Hubbard model on a square lattice falls 
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in the first category (i), since its sign problem is entirely due to the Fermi 
statistics of the electrons [23]. 

2.2.2. Importance sampling and Fixed-Node approximation 
As clearly explained in at least one of the other chapters of this book [17], 
the use of importance sampling is of crucial importance to improve the sam- 
pling efficiency; a trial function may also be the simplest way to incorpo- 
rate the Fermi statistics in a basis of "labeled fermion" lattice configurations 
[13, 23]. So in actual fermion calculations a good trial function ^>t (typically 
available from preliminary Variational Monte Carlo calculations) is always 
adopted. In the power method, importance sampling amounts to replac- 
ing g by g, whose matrix elements are g X 'x = ^T(X')(X'\g\X)^ 1 (X), 
and then to consistently deal with g rather than g [17]. The sign problem 
amounts then to the occurrence of some unwanted gx'X < 0, and, apart 
from transient estimates which are rather hazardous for large systems, the 
only cure to this problem seems to be the fixed-node approximation re- 
cently proposed by An, Bemmel, Ceperley, Haaf, Leeuwen, and Saarloos 
[13, 24]. In this approximation the true hamiltonian Ti is replaced by an 
effective hamiltonian 7i e f[, defined as follows: 

{X'\H eS \X) = ifg X >x<0 

= (X'\H\X) (otherwise) (1) 

for the off-diagonal terms X ^ X', so that hops from the state \X) towards 
"sign-flipping states" \X') (i.e. such that g < 0) are now forbidden, and, 
for the diagonal terms, 

(x\n cS \x) = (x\n\x)+Y,(x'\n\x)^^, (2) 

the summation being over all the "sign-flipping states" \X') such that 
(X'\H\X) but g X 'x < 0: whenever in the "neighborhood" of the 
state \X) there are sign-flipping states \X'), appropriate repulsive terms 
are added to the old diagonal element (X\H\X). Such a diagonal repulsive 
term has the intuitive effect of repelling the random walk away from the 
regions of the Hilbert space where g would tend to revert its sign, and 
also the less intuitive, but more important effect of producing, in analogy 
with the continuum case [25, 26], an upper bound for the true ground-state 
energy and a variational approximation [13]. This formulation, contained 
in Ref.[13], extends the original prescription of Ref. [24] to include all the 
possible sources of a sign-flip in g (Subs. 2.2.1) [27]. However for our sim- 
ple Hubbard hamiltonian (Eq. 3) the original and new prescription are 
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coincident, because in this case the only source of sign flips are the Fermi 
nodes [23]. Although the original papers by An, Bemmel, Ceperley, Haaf, 
Leeuwen, and Saarloos [13, 24] are only 3-4 years old, their "fixed-node 
Monte Carlo for lattice fermions" has been already used by other authors, 
like Boninsegni for frustrated Heisenberg systems [22], or Koch et al. for 
orbitally-degenerate Hubbard models [4, 28]. We (Cosentini et al. [11]) have 
chosen it for the extensive study of the 2D Hubbard model presented in the 
coming Sections 3,4, because it allows a reliable investigation of (previously 
unfeasible) very large lattice-fermion systems. 

3. 2D Hubbard model and phase separation 

3.1. MOTIVATIONS 

3.1.1. High-T c , phase separation, ICDW 

Strongly correlated electrons and holes are expected to play a key role in 
the high-T c superconductors. Their possible instability towards phase sep- 
aration, initially believed to inhibit superconductivity, is attracting a lot of 
interest since a few different authors [29, 30, 31] have pointed out that such 
a tendency may in fact be intimately related to the high-T c superconductiv- 
ity. Long-range repulsive interactions may turn the phase-separation insta- 
bility into an incommensurate charge-density-wave (ICDW) instability, and 
the very existence of a quantum critical point associated with it may be a 
crucial ingredient of the superconducting transition [32]. Phase-separation 
and/or ICDW instabilities are related to a substantial reduction of the ki- 
netic energy, which otherwise tends to stabilize uniformly distributed states; 
such a reduction is typical of strongly correlated electrons, both in real and 
model systems. 

3.1.2. What is phase separation? 

Phase separation is a thermodynamic instability associated with the vio- 
lation of the stability condition x -1 = d 2 £/dn 2 > 0, which requires the 
energy density 6 of an infinite electronic system to be a convex function 
of the electron density n. If such a violation occurs in the density range 
n\ < n < ri2, then the system will separate into two subsystems with elec- 
tron densities n\ and n<i. For the two-dimensional t—J and Hubbard models 
a phase separation, if any, is expected to occur in a density range close to 
half filling (n ~ 1), and to yield a hole-rich phase with electron density 
ni < 1 and a hole-free phase with electron density ri2 = 1 [33]. In view of 
our further developments, the same physics can be conveniently rephrased 
in terms of the hole density x = 1 — n and energy density e/j(x) = £(n(x)), 
as done in the upper right panel of Fig. 1. 
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3.1.3. Maxwell's construction 

In a truly infinite system a phase separation of the type just discussed 
would be associated to a vanishing inverse compressiblity x~ l i n the whole 
density range n\ < n < ri2\ in a finite system \ l may even become nega- 
tive, because of surface effects; the same applies of course to the the second 
derivative of the hole energy density e/ l (x) with respect to the hole density 
x. So for finite systems it's preferable to pinpoint the phase separation us- 
ing Maxwell's construction (Fig. 1, upper left panel). It has been shown in 
Ref . [33] that Maxwell's construction is equivalent to study, as a function of 
the hole density x = l— n, the quantity e(x) = [e/j(x) — en]/x, i.e. the energy 
per hole (x) measured relative to its value at half filling en = ey l {x = 0) 
(Fig. 1, lower panels). For an infinite system, if the inverse compressibility 
X^ 1 vanishes between some critical density x c and half filling x = (up- 
per right), then the function e(x) is a constant for < x < x c , and the 
fingerprint of a phase separation is thus a horizontal plot of e(x) below x c 
(lower right). For a finite system, instead, the inverse compressibility x -1 
may become negative (upper left) and the fingerprint of a phase separa- 
tion is a minimum of e(x) at x = x c [33] (lower left). On finite lattices, 
where only a discrete set of densities is available and a statistical error bar 
usually accompanies the estimated QMC energies, the plot of the function 
e{x) (lower panels of Fig. 1) is supposed to allow an easy visual judgement 
of the presence or absence of such an instability. The attempt of telling 
whether near half filling a few electron or hole energies are better interpo- 
lated by a straight or by a somewhat curved line (upper panels of Fig. 1) is 
definitely more difficult than looking at e(x). However even e(x) can give a 
reliable criterion only for medium-large finite systems; really small systems 
(for which most numerical results have been up to now available) can at- 
tain so few and coarse densities, and suffer from so large finite-size errors, 
that their predictions of the relevant trends remain largely inconclusive, 
no matter which energy function we look at. More remarks on e(x) vs. 
phase separation will be presented later on, when examining the numerical 
results. 

3.1.4. Limited experimental evidence 

Phase separation has been experimentally observed in La2Cu04 + 5[34, 35], 
where the oxygen ions can move: in the doping interval 0.01 < 5 < 0.06 the 
compound separates into a nearly stoichiometric antiferromagnetic phase 
and a superconducting oxygen-rich phase. In generic compounds, where 
charged ions cannot move, the possibility of a macroscopic phase separa- 
tion is spoiled by the long-range Coulomb repulsion, and should lead to an 
ICDW instability [36]; here, however, the identification of charge inhomo- 
geneities with spoiled phase separation is less obvious [37]. 
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Figure 1. Sketch of a hypothetic system of holes which, when its density falls between 
zero and a critical value x c , tends to phase-separate between a hole- free phase (zero 
density) and another hole-rich phase (density equal to x c ). Upper panel: behavior of the 
hole energy density en vs. the hole density x in the case of a finite (left) and an infinite 
(right) system. For the finite system (left), between x — and x = x c , the actual energy 
of the phase-separated system (dashed straight line, Maxwell's construction) is lower 
than the energy of the homogeneous system; for the infinite system, in the same density 
range, the energy density is a straight line, and the inverse compressibility vanishes. 
Lower panel: behavior of the function e(x) — [eh{x) — en]/ x for a finite (left) and an 
infinite (right) system. As shown by Emery et at, in a finite system the study of e(x) is 
completely equivalent to Maxwell's construction (see text). 



3.1.5. Theoretical controversy 

On the theoretical side, evidence for phase separation has been suggested 
for various models of strongly correlated electrons, such as the t—J model 
[33], the three-band Hubbard model, the Hubbard- Holstein model and the 
Kondo model (see e.g. Ref. [32] and references therein). Despite intensive 
studies, even for simple models there is no general agreement on the phase 
separation boundary: for the very popular t — J model, phase separation 
is fully established only at large J, but at small J (which unfortunately 
happens to be the physically relevant case) theoretical and numerical results 
are quite controversial. Emery et aVs [33] theory that phase separation 
occurs at any value of J in the t—J model is confirmed by a recent numerical 
study by Hellberg and Manousakis [38], but is in contrast with Dagotto et 
al.'s [31] exact numerical results on small clusters, suggesting no tendency 
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toward phase separation for both the Hubbard model and the t — J model 
below a critical value J < J c ~ t, and with Shih et al.'s [39] numerical 
results. We also mention the recent suggestion by Gang Su [40], according 
to which the Hubbard model does not show phase separation for any value 
of U/t at any temperature. 

3.1.6. A challenging application of fixed-node quantum Monte Carlo 
Given the above theoretical and experimental situation, the availability of 
the fixed-node quantum Monte Carlo, described in the above Subs. 2.1, 
has recently provided us (Cosentini et al. [11]) with a strong motivation to 
further investigate the Hubbard model. Whether the 2D Hubbard hamil- 
tonian, a prototype for interacting electrons with no long-range repulsion, 
shows any instability towards phase separation, is by itself an interesting 
open question. Besides the intrinsic interest in the model's behavior, a nu- 
merical investigation of this question may also shed some indirect light on 
two related issues: (i) whether phase separation is likely to occur in the 
t — J model at small J (the "physical region", where the ground-state en- 
ergy difference with the Hubbard model should become negligible); and 
(ii) whether the measured phase separation of some real high-T c supercon- 
ductors (Subs. 3.1.4) is such a fundamental and elementary fact as to be 
explained by the simplest one-band Hubbard hamiltonian. 

3.2. NUMERICAL STRATEGY AND TESTS 

3.2.1. Hamiltonian 

To obtain numerical estimates of functions like those sketched in Fig. 1 we 
need to to evaluate the ground-state energy of the Hubbard hamiltonian 

W = "* E ( C *W + h - c -) + U ]T n t ^n tl (3) 

for many large finite lattices and many different electron densities. As al- 
ready mentioned (and easily double-checked by the reader), for fermions 
this hamiltonian H and the corresponding importance-sampled Q satisfies 
the requirement of sparseness (Subs. 2.1.1 and Ref. [18]) but not the one on 
sign (Subs. 2.2 and Ref. [23]). The choice of a good variational wavefunction 
is thus needed not only to improve the sampling efficiency, but also for the 
power method to become feasible through the fixed-node approximation 
(Subs. 2.2.2). 

3.2.2. Choice of variational wavefunction 

The variational wavefunctions we use to guide the random walks and to 
fix the nodes are the product of a Gutzwiller factor and two Slater deter- 
minants of single-particle, mean-field wavefunctions (with uniform electron 
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density n and uniform staggered magnetization density m) for an equal 
number of up- and down-spin electrons [41]. Preliminary estimates of the 
optimal Gutzwiller parameter < g < 1 and mean-field wavefunctions 
(parametrized by the staggered magnetization < m < 1) were obtained, 
for each U/t and density, by variational Monte Carlo (VMC) runs; they are 
shown in Fig. 2 as a function of the electron density n for three choices of the 
coupling strength U (in units of the hopping parameter t). As U ft increases, 
the optimal g decreases, reflecting a stronger depression of doubly-occupied 
configurations. As far as the density dependence is concerned, the staggered 
magnetization m is close to zero for all U ft and low-medium electron den- 
sities, but sharply raises above n ~ 0.75, thus signaling, at the variational 
level, the system's tendency towards an antiferromagnetic ordering as half 
filling approaches. The optimal Gutzwiller factor steadily decreases as a 
function of the electron density n, except for very high electron densities 
(n ~ 0.9 and above), where it apparently goes up again, as also noticed 
elsewhere [42] . In this narrow density range, however, the energy minimum 
appears to be flat: slightly different g values give essentially the same VMC 
energy, and the corresponding trial wavefunctions also yield the same en- 
ergy at the fixed-node level. 

3.2.3. Comparison with exact results 

A few representative variational (VMC) and fixed-node Monte Carlo (FNMC) 
energies are shown in Table 1 for the 4x4 Hubbbard lattice, for which exact 
results [43] are available. As expected, the VMC energy is always above the 
fixed-node energy, which in turn, for these coupling strengths, is slightly 
(~ 3%) above the exact energy. For comparison we show the Constrained- 
Path Monte Carlo (CPMC) energies of Zhang et al. [44, 45], which also 
include a larger 16x16 lattice (last row). Especially at large U/t our results 
appear to be of comparable quality to theirs. As far as the 4x4 results are 
concerned, we notice that for N e = 10, which corresponds to a closed-shell 
configuration, both fixed-node and CPMC are much closer to the exact en- 
ergy than for N e = 14, which corresponds to an open-shell configuration. 
This could be a serious problem when numerically studying the behavior 
of the energy as a function of the density; the results presented here sug- 
gest that, for lattices larger than 12 x 12, the shell effects become almost 
irrelevant [47]. 

3.2.4. How to vary the density? 

To study the energy as a function of the electronic density we have first 
tried out the less usual way of varying the density suggested in Ref. [38] to 
"avoid spurious Fermi-surface shape effects" : keep the number of electrons 
N e fixed while the size of the underlying lattice is varied. But we discovered 
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Figure 2. A collection of optimal Gutzwiller factors < g < 1 (upper panel) and 
staggered magnetizations < m < 1 (lower panel) for many different electron densities 
and three values of the coupling strength U/t: 10 (triangles), 20 (diamonds), and 40 
(dots). Most results correspond to 16x16 lattices, but a subset corresponds to 12x12 
and to other large but different lattices. When the size is large enough, the optimal g,m 
appear thus to be almost independent of the size. 



size 


N e 


n 


U/t 


VMC 


FNMC 


CPMC 


EXACT 


4x4 


10 


0.625 


4 


-1.211(2) 


-1.220(2) 


-1.2238(6) 


-1.2238 


4x4 


10 


0.625 


8 


-1.066(2) 


-1.086(2) 


-1.0925(7) 


-1.0944 


4x4 


14 


0.875 


8 


-0.681(2) 


-0.720(2) 


-0.728(3) 


-0.742 


4x4 


14 


0.875 


12 


-0.546(2) 


-0.603(2) 


-0.606(5) 


-0.628 


16x16 


202 


0.789 


4 


-1.096(2) 


-1.107(5) 


-1.1193(3) 





TABLE 1. Ground-state energy per site (in units of the hopping parameter t) for a 
4x4 Hubbard lattice and various values of U/t. N e is the number of electrons and n 
is the corresponding average density. VMC: variational Monte Carlo, Ref. [11]; FNMC: 
Fixed-Node Green's function Monte Carlo, Ref. [11]; CPMC: Constrained-Path Monte 
Carlo, Ref. [44]; EXACT: exact diagonalization results, Ref. [43]. (see text). 
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that by this prescription, if the number of electrons is really small (e.g. 
N e = 16), then artificial changes in the convexity of the curve may occur. 
If, instead, the system is large enough (e.g. 12 x 12 lattices or larger), then 
it doesn't matter how the density is varied, as this prescription and the 
usual prescription give the same results. So for our systematic study (many 
densities and three U/t values) we finally adopted a large 16x16 lattice 
(N s = 256 sites), and varied the number of electrons N e to yield electronic 
densities n = N e /N s ranging from empty n = to half filling n = l. 

3.3. RESULTS 

3.3.1. Energy vs. density for large lattices 

In Figs. 3-6 we show the electronic ground-state energy per site, obtained by 
FNMC runs as a function of the electron density [46] . Energies are in units 
of the hopping parameter t throughout this paper; the statistical errors 
are smaller than the marker size, and thus are not visible. The calculated 
points are shown as full markers for closed shells, and as empty markers for 
open shells. From the comparison of Fig. 3 and Table 1 it appears that the 
open-shell error, which we found to be significant for a small 4x4 lattice, 
becomes of the order of the statistical error (and thus negligible [47]) for 
our large lattices [48]. 

3.3.2. Accuracy 

At all densities our three sets of data for U/t = 10 (Fig. 3), 20 (Fig. 4), and 
40 (Fig. 5) are bracketed, as seen in the cumulative Fig. 6, by the noninter- 
acting unpolarized energy and the fully spin-polarized energy (both dashed 
in Fig. 6), and display a smooth and reasonable behavior. To evaluate the 
absolute accuracy of our results, we can rely on two exact limits: the low- 
density (n~0) regime, where we expect £ = — 4n, and the half filled case 
(n = l), for which the strong-coupling expansion provides the correct large- 
U/t behavior: to leading order in t/U, the model maps onto an Heisenberg 
model, whose ground-state energy has been evaluated with great accuracy 
[12, 49]. We can also consider the next correction term 34.6t 4 /[/ 3 [50]. At 
low density our results are essentially exact; at half filling our error is small 
(~ 3%) for U/t = 10 but (as already noticed in Table 1) it tends to grow 
with U/t: for U/t = 20 it's ~9%, and for U/t = 40 it's ~ 11%. We have made 
sure (see markers other than dots in Figs. 3,5 and Refs. [46, 48]) that such 
an energy discrepancy is not due to shape, open-shell, finite-size [47], and 
boundary condition effects; as far as systematic errors are concerned, we 
are thus left with the fixed- node approximation: as U/t grows, more flexible 
trial wavefunctions may be required to obtain more accurate energies [51]. 
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Figure 3. Ground-state energy per site (in units of the hopping parameter t) as a 
function of the electronic density, for a 2D Hubbard lattice of N a = 16 x 16 = 256 sites 
with U/t — 10. Errors are smaller than the dot size. Full dots correspond to closed shells 
and empty dots correspond to open shells. The three crosses correspond to closed-shell 
densities of a ll\/2 x llv^ lattice, and are shown for comparison (see text and Refs.[46, 
48]). 
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Figure 4- Ground-state energy per site (in units of the hopping parameter t) as a 
function of the electronic density, for a 2D Hubbard lattice of N„ — 16 x 16 = 256 sites 
with U/t — 20. Errors are smaller than the dot size. Full dots correspond to closed shells 
and empty dots correspond to open shells. 
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Figure 5. Ground-state energy per site (in units of the hopping parameter t) as a 
function of the electronic density, for a 2D Hubbard lattice of N s — 16 x 16 = 256 sites 
with U/t — 40. Errors are smaller than the dot size. Full dots correspond to closed shells 
and empty dots correspond to open shells. Full and empty triangles, corresponding to 
closed and open shells of a smaller 12 x 12 lattice, are shown for comparison, and give a 
measure of the size effects (see text and Ref. [47]). 

3.3.3. Phase separation 

Keeping in mind the virtues and limitations of our numerical study, we 
can now turn to the question of phase separation in the Hubbard model. 
In some sense the function e(x) works like a magnifying lens of the phase 
separation (see Fig. 1). It should be stressed that in a consistent definition 
of e{x) the half-filling energy en must be obtained as e^ix = 0) from the 
same calculation as any other e^(x 7^ 0) (here, from the FNMC at half 
filling). If that's not the case (for example, if the Heisenberg value is used 
instead), then e(x) will spuriously diverge near x = 0, with a good chance 
of artificially creating, rather than magnifying, the occurrence of phase 
separation. In Figs. 7-9 we find plots of e(x) for U/t = 10, 20, and 40; these 
values, as well as the associated error bars, are directly obtained from those 
of Figs. 3-5 (i.e. from the original FNMC energies and tiny error bars). 
Despite the error bars, a common trend is evident for all the calculated 
coupling strenghts: e{x) has a positive slope for large hole densities, far 
from half-filling, but then it clearly changes slope below some small critical 
density x c . Such a minimum in e(x) implies that, at least for the FNMC 
effective hamiltonian determined by our choice of wavefunction and cell 
size [47], phase separation occurs below x = x c . Although a finer grid of 
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Figure 6. This picture summarizes the results of the previous Figs. 3-5, showing the 
ground-state energy per site (in units of the hopping parameter t) as a function of the 
electronic density, for a 2D Hubbard lattice of N s = 16x16 = 256 sites with U/t = 10 
(lower), 20 (middle), and 40 (upper data). Errors are smaller than the dot size and 
symbols are the same as in Figs. 3-5. The dashed curves correspond to two known results 
for the infinite lattice at U — 0: the fully spin-polarized case (upper curve), whose total 
energy per site is symmetric with respect to quarter filling, and the unpolarized case 
(lower curve), whose total energy per site is symmetric with respect to half filling, (see 
text). 

hole densities would be required to locate with high precision the critical 
density function of U/t, we already see that x c decreases as U/t 

is increased; this qualitatively agrees with the original predictions [33] and 
with some previous calculations on the t—J model at corresponding values 
of J = At 2 /U [38]. 

4. Summary, open problems, conclusions and perspectives 

4.1. SUMMARY 

In summary, the extensive fixed-node Monte Carlo simulations of the Hub- 
bard model for 16x16 two-dimensional lattices by Cosentini et al. [11], just 
recalled in the previous Section, suggest the presence of a phase separation 
for U 3> t. If confirmed, this result would imply that the t—J model is also 
likely to show a phase separation in the physically relevant regime J < 0.4, 
and that even a single-band Hubbard model is sufficient to reproduce this 
physical tendency of some high-T c superconductors (Subs. 3.1). Before these 
statements are confirmed beyond any conceivable doubt, further tests may 
be required because of the following two problems. 
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Figure 7. Plot of e(x) vs. cc for U/t = 10. Full dots correspond to closed shells and empty 
dots correspond to open shells of a 16x16 lattice. Crosses (corresponding to a llv2x llv2 
lattice) are shown for comparison. Obviously at small x the error bar associated to e(x), 
Ae(x) = [Aeh(x) + Aeh(x = 0)]/x, becomes significant even if the statistical FNMC error 
Aeh{x) is tiny. 
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Figure 8. Plot of e(x) vs. x for U/t = 20. Full dots correspond to closed shells and 
empty dots correspond to open shells of a 16 x 16 lattice. Other comments like in the 
previous figure. 

4.2. OPEN PROBLEMS 

4.2.1. Fixed-node error at large U/t 

First of all, the fixed-node approximation is a variational approximation. 
Especially for the two larger U/t values, the energy discrepancy at half 
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Figure 9. Plot of e(x) vs. x (see text) for U/t = 40. Full dots correspond to closed 
shells and empty dots correspond to open shells of a 16 x 16 lattice. Full and open trian- 
gles correspond to closed and open shells of a smaller 12 x 12 lattice and are shown for 
comparison (see text). Other comments like in the previous figure. 

filling (as discussed in Subs. 3.3.2) is unfortunately significant. Not knowing 
in advance the density dependence of such an energy discrepancy, even the 
conclusions on phase separation, based on the function e(x), could be at 
risk. Richer variational wavefunctions (yielding different nodal topologies) 
[51], or possibly the stochastic reconfiguration scheme recently proposed by 
Sorella [53], must be employed to go beyond our simple Gutzwiller-Slater 
nodes, and to settle this point completely. 

4.2.2. Finite-size and "phase separation" at U = 

A second motivation for further tests is illustrated in Fig. 10. Here we show 
the energy density of a finite 16x16 lattice with U/t = (noninteracting 
electrons). Close to half filling, below x c = 0.125, the function e(x) becomes 
completely flat. Such a behavior, already pointed out by Lin a few years 
ago [54] but surprisingly forgotten in most subsequent numerical studies of 
small Hubbard lattices, is an artifact related to the high degeneracy of the 
Fermi level, and disappears in a truly infinite lattice (x c vanishes as 2/L 
when L, the side of a square LxL lattice, goes to infinity). Comparing Fig. 10 
with previous Figs. 7-9, one may even fear that the artificial flattening at 
U = has something to do with the behavior found for e(x) at U > 0. This 
observation suggests caution even if, as stated in various parts of Subs. 3.3, 
many severe tests of open-shell and finite-size effects have already been 
passed. The almost perfect convergence of energies £(n), shown in Figs. 3- 
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Figure 10. This figure contains plot of e(x) vs. x for Z7/t = (noninteracting electrons, 
analytic results) for a 16 x 16 lattice. As a marker we use here vertical bars; we display 
data for all the 128 densities which can be obtained at this lattice size with an even 
number of electrons. Note that below x = 2/16 e(x) is completely flat. 



5, and even the convergence (within the appropriate statistical error) of the 
e(x) data of Figs. 7-9, could, incredibly, not be enough to rule out a very 
slow 2/L drift of the critical density towards zero. Perturbation theory or 
further numerical work at small U/t, as well as a thorough finite-size scaling 
study based on yet larger lattices may be required to convincingly rule out 
the possibility of such a spurious effect. 

4.3. CONCLUSIONS AND PERSPECTIVES 

Neither previous, quite crude numerical studies, nor the recent (and defi- 
nitely much more accurate) work by Cosentini et al, recalled in this chapter, 
seem to be sufficient for a 100% safe and general claim on phase separation. 
This is by itself information of the greatest importance to whoever is active 
either in the theory or in the numerical simulation of strongly correlated 
two-dimensional electron systems: phase separation is a very delicate effect, 
which can be masked by otherwise negligible variational or finite-size errors. 
The simplicity and ease of the lattice fixed-node Monte Carlo [24, 13], the 
efficiency of the new Calandra-Sorella reconfiguration scheme [12], and the 
experience gained in the preliminary study recalled here [11] provide a good 
starting point and very useful guidelines for a future, conclusive numerical 
study of this challenging problem. 
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